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AN ANALYSIS OF THE LI SCHEME FOR THE SUBDIFFUSION 
EQUATION WITH NONSMOOTH DATA 


BANGTI JIN, RAYTCHO LAZAROV, AND ZHI ZHOU 


Abstract. The subdiffusion equation with a Caputo fractional derivative of order 
OL G (0, 1) in time arises in a wide variety of practical applications, and it is often 
adopted to model anomalous subdiffusion processes in heterogeneous media. The 
LI scheme is one of the most popular and successful numerical methods for dis¬ 
cretizing the Caputo fractional derivative in time. The scheme was analyzed earlier 
independently by Lin and Xu (2007) and Sun and Wu (2006), and an 0(r^““) 
convergence rate was established, under the assumption that the solution is twice 
continuously differentiable in time. However, in view of the smoothing property of 
the subdiffusion equation, this regularity condition is restrictive, since it does not 
hold even for the homogeneous problem with a smooth initial data. In this work, 
we revisit the error analysis of the scheme, and establish an 0(t) convergence rate 
for both smooth and nonsmooth initial data. The analysis is valid for more gen¬ 
eral sectorial operators. In particular, the LI scheme is applied to one-dimensional 
space-time fractional diffusion equations, which involves also a Riemann-Liouville 
derivative of order ^ G (3/2,2) in space, and error estimates are provided for the 
fully discrete scheme. Numerical experiments are provided to verify the sharpness 
of the error estimates, and robustness of the scheme with respect to data regularity. 
Keywords: fractional diffusion, LI scheme, error estimates, space-time fractional 
diffusion 


1. Introduction 

We consider the model initial-boundary value problem for the following fractional order 
parabolic differential equation for u{x,t): 

dfu — Au = /, in fl T > t > 0, 

(1.1) u = 0, on 90 T>t>0, 

u(0) = V, in O, 

where O is a bounded convex polygonal domain in R'^ (d = 1, 2, 3) with a boundary 90, 
u is a given function on O, and T > 0 is a fixed value. Here d^u (0 < a < 1) denotes 
the left-sided Caputo fractional derivative of order a with respect to t and it is defined by 
(see, e.g. [13, pp.91]) 

(1.2) 9t“u(t) = ——1— -f {t-s)~°‘u{s)ds, 

r(i - a) Jo 

where r(-) is Euler’s Gamma function defined by r(a;) = 

The model (1.1) is known to capture well the dynamics of subdiffusion processes, in 
which the mean square variance grows at a rate slower than that in a Gaussian process [4] , 
and has found a number of applications. For example, subdiffusion has been successfully 
used to describe thermal diffusion in media with fractal geometry [19], highly heteroge¬ 
neous aquifer [2[ and underground environmental problem [7[. At a microscopic level, the 
particle motion can be described by a continuous time random walk, in which the waiting 
time of the particle motion follows a heavy tailed distribution, as opposed to a Gauss¬ 
ian process, which is characteristic of the normal diffusion equation. The macroscopic 
counterpart is a diffusion equation with a Caputo fractional derivative in time, i.e., (1.1). 
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The derivation and study of accurate numerical methods for the model (1-1) with 
provable (possibly optimal-order) error bounds have attracted considerable interest in 
recent years, and a number of efficient numerical schemes, notably the finite difference 
method, have been developed. There are two predominant approximations: the LI type 
approximation and Grunwald-Letnikov type approximation. The former is essentially of 
finite difference nature, whereas the latter is often based on convolntion quadrature. We 
refer interested readers to [11, Section 1 and Table 1] for an updated overview of these 
numerical methods and relevant convergence rate results. 

In this paper, for reasons to be explained below, we revisit the LI scheme, which was 
independently developed and analyzed in [22] and [16]. To this end, we divide the interval 
[0, T] into a uniform grid with a time step size t = T/N, N G N, so that 0 = to < ti < 

... < In = T, and = nr, n = 0,..., N. The LI approximation of the Caputo fractional 
derivative dtu{x,tn) is given by [16, Section 3] 


dtu{x,t„) = 


r(i- 


1 y [ 


du{x, s) 


1=0 •' 
n — 1 


ds 


(tn — S) ds 


(1.3) 


r(l-a) 




{x,tj+i) - u{x,tj) r^+ 


1=0 


/■L+1 

Jti 


{tn — S) °‘ds 




u{x, tn-j) - u{x, tn-j-l) 


1 = 0 


= r °‘[bou{x,t„) - bn-iu{x,to) + ^(&i - bj-i)u{x,tn-j)] =■■ Li{u). 

1=1 


where the weights bj are given by 

bj = iU + !)'■“ - /■“)/r(2 - o), j = 0,1,. .., iV - 1. 


It was shown in [16, equation (3.3)] (see also [22, Lemma 4.1]) that the local truncation 
error of the LI approximation is bounded by cr^““ for some constant c depending only 
on u, provided that the solution u is twice continuously differentiable. Further, upon dis¬ 
cretizing the spatial derivative(s) using the finite difference method, finite element method 
or spectral method, we arrive at a fully discrete scheme. Since its first appearance, the LI 
scheme has been extensively used in practice, and currently it is one of the most popular 
and successful numerical methods for solving the time fractional diffusion equation (1.1), 
including the case of nonsmooth data arising in inverse problems (see, e.g., [12]). 

Note that the fractional diffusion operator has only very limited smoothing property, 
especially for t close to zero. For example, for the homogeneous equation with an initial 
data V € L^{Q), we have the following stability estimate [20, Theorem 2.1] 


< ct “|]u]]i 2 (n). 

That is, the a-th order Caputo derivative is already unbounded, not to mention the high- 
order derivatives. In view of the limited regularity of the solution of (1.1), the high 
regularity required in the convergence analysis in these useful works is restrictive. 

Hence the C^-regularity assumption generally does not hold for problem (1.1), and 
the case of nonsmooth data is not covered by the existing error analysis. Our numerical 
experiments indicate that the 0(r^““) convergence rate actually does not hold even for 
smooth initial data v. This is clearly seen from the numerical results in Table 1, where 
rate denotes the empirical convergence rate, and the numbers in the bracket denote the 
theoretical rates based on the local truncation error. These computational results show 
that the scheme is only first-order accurate even when the initial data v is smooth. This 
observation necessitates revisiting the convergence analysis of the LI scheme, especially 
for the case of nonsmooth problem data. 
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The goal of this paper is to fill the gap between the existing convergence theory and the 
numerical experiments, namely, establishing of optimal error bounds that are expressed di¬ 
rectly in terms of the regularity of the problem data. For the standard parabolic equation, 
such bounds are well known (see, e.g. [23, Chapter 7]). To the best of our knowledge, 
there is no such error analysis of the LI scheme in the case of nonsmooth data. Such 
theory and numerical experiments are provided in this work. 

Table 1. The L^-norm of the error ||17^ ~i^(in)||L 2 (n) for problem (1.1) 
with / = 0 and (a) v = X(o,i/ 2 ) and (b) v = *(1 — x), at t = 0.1, 
computed by the fully discrete scheme (2.8) with h = 2~ . 


a 

N 

10 

20 

40 

80 

160 

320 

rate 

0.1 

(a) 

3.30e-4 

1.62e-4 

8.02e-5 

3.99e-5 

1.99e-5 

9.93e-6 

« 1.01 (1.90) 


(b) 

4.97e-4 

2.44e-4 

1.21e-4 

6.00e-5 

2.99e-5 

1.49e-5 

« 1.00 (1.90) 

0.5 

(a) 

3.04e-3 

1.44e-3 

6.96e-4 

3.41e-4 

1.68e-4 

8.35e-5 

« 1.03 (1.50) 


(b) 

4.61e-3 

2.18e-3 

1.05e-3 

5.16e-4 

2.54e-4 

1.34e-5 

« 1.02 (1.50) 

0.9 

(a) 

1.31e-2 

6.46e-3 

3.17e-3 

1.56e-3 

7.68e-4 

3.78e-4 

« 1.02 (1.10) 


(b) 

1.94e-2 

9.67e-3 

4.78e-3 

2.36e-3 

1.16e-3 

5.75e-4 

« 1.02 (1.10) 


In Theorem 3.3, we present an optimal 0(r) convergence rate for the fully discrete 
scheme based on the LI scheme (1.3) in time and the Galerkin finite element method in 
space for both smooth and nonsmooth data, i.e., v € L^(fl) and Av € L^{Q) {A = — A with 
a homogeneous Dirichlet boundary condition), respectively. For example, for v £ L^{Q) 
and ( 7 ° = PhV, for the fully discrete solution ? 7 ^, there holds 

||u(U) - [/;;‘||i2(n) < c{Tt~^ + h'^t~°‘)\\v\\L2(^ny 

Surprisingly, for both v € L^{Q,) and Av G L^(fl), the error estimate deteriorates as time 
t approaches zero, but for any fixed time tn > 0, it can achieve a first-order convergence. 
Extensive numerical experiments confirm the optimality of the convergence rates. Our 
estimates are derived using the techniques developed by [17] for convolution quadrature 
and in the interesting recent work of [18] on a piecewise constant discontinuous Galerkin 
method. The proof essentially boils down to some delicate estimates of the kernel function, 
which involves the poly logarithmic function. Finally, we note that our results are appli¬ 
cable to more general sectorial operators, including the very interesting case of the space 
time fractional differential problem involving a Riemann-Liouville derivative in space [9]. 

The rest of the paper is organized as follows. In Section 2, we recall preliminaries on the 
fully discrete scheme, and derive the solution representation for the semidiscrete and fully 
discrete schemes, which play an important role in the error analysis. The full technical 
details of the convergence analysis are presented in Section 3. In Section 4 we consider the 
adaptation of the scheme to a one-dimensional time-space fractional differential equation, 
and derive optimal convergence rate. It appears to be the first error estimates expressed 
directly in terms of the data regularity for such equation with nonsmooth data. Numerical 
results are presented in Section 5 to confirm the convergence theory and the robustness of 
the scheme. Throughout, the notation c, with or without a subscript, denotes a generic 
constant, which may differ at different occurrences, but it is always independent of the 
mesh size h and the time step size r. 

2. Preliminary 

In this part, we give the semidiscrete and fully discrete schemes, based on a standard 
Galerkin method in space and the LI approximation in time. 
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2.1. Semidiscrete scheme. Since the solution u : (0,T] —>■ L^{yi) can be analytically 
extended to the sector {z G C; z ^ 0, |arg 2 :| < 7r/2} [20, Theorem 2.1], when / = 0, we 
may apply the Laplace transform to equation (1.1) to deduce 

(2.1) z°‘u{z) + Au{z) = z°‘~^v, 

with the operator A = —A with a homogeneous Dirichlet boundary condition. Hence the 
solution u{t) can be represented by 

( 2 . 2 ) = e^\z°‘I+ A)-^z°‘-\dz, 

where the contour Fa,a is given by 

re,« = {z G C : \z\ = S,\ argzj < 9} yj {z G C ■. z = , p > 5}. 

Throughout, we choose the angle 9 G (7r/2,7r). Then G Eg/ with 9' — a9 < -k for all 
z G Tie '■= {z G C ) \ argzj < 9}. Then there exists a constant c which depends only on 9 
and a such that 

(2.3) \\{z°‘I + A)-^\\<cz-°‘, VzG^e. 

Now we introduce the spatial semidiscrete scheme based on the Galerkin hnite element 
method. Let Th he a. shape regular and quasi-uniform triangulation of the domain H 
into d-simplexes, denoted by T. Then over the triangulation Th we define a continuous 
piecewise linear finite element space Xh by 

Xh = {vh G Hq{Q,) : Vh\T is a linear function, VT GTh} ■ 

On the space Xh, we define the (f2)-orthogonal projection Ph : L‘^{Q) —>■ Xh and the 
Ritz projection Rh : Hq (Q,) —>■ Xh , respectively, by 

{PhV>,x) = i‘P,x) VyeAh, 

{VRh<fi,Xx) = {'Vip,Xx) VxeXh, 

where (•, •) denotes the Z/^(H)-inner product. Then the semidiscrete Galerkin scheme for 
problem (1.1) reads: Hnd Uh{t) G Xh such that 

(2.4) id:uh,x) + (Vm^, Vx) = if, X) yxeXh, 

with Uh(0) = Vh G Xh- Upon introducing the discrete Laplacian Xh : Xh —>■ Xh defined 
by 

-(Aj.v5,x) = (Vv?, Vx) Vv3, xeXh, 
the spatial semidiscrete scheme (2.4) can be rewritten into 

(2.5) dtUhit) + AhUhit) = fhit), t> Q 

with UhiO) = Vh G: Xh, fh = Phf and Ah = —Ah- Like before, the solution Uh to (2.5) 
with fh = 0 can be represented by 

(2.6) Wh(i) = 7^ / e‘'\z°‘+ Ah)~^z°‘~^Vhdz- 

Further, for later analysis, we let Wh = Uh — Vh- Then Wh satisfies the problem: 

dtWh + AhWh = -AhVh, 
with Wh{0) = 0. The Laplace transform gives 

z°‘wh{z) + AhWhiz) = -z~^AhVh- 
Hence, Wh{z) = K\{z)vh, with 

Ki{z) = -z-\z’^I + Ah)-^Ah, 
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and the desired representation for Wh{t) follows from the inverse Laplace transform 

(2.7) Whit) =-^ [ e^*Ki{z)vhdz. 

The semidiscrete solution Uh satisfies the following estimates from [10]. The log factor 
in the estimates in [10, Section 3] can be removed using the operator trick in [3, Section 
3]. 

Theorem 2.1. Let u and Uh be the solutions of problems (1.1) and (2.5) with / = 0 and 
fh = 0, respectively. Then the following error estimates hold. 

(a) If Av € L^{Q,) and Vh = RhV, then 

\\uit) - Uh(t)||i 2 (n) < c/i^||7ln||i2(n). 

(b) If V & L^{Ti) and Vh = PhV, then 

IWitn) - Wh(t)|L2(n) < 

2.2. Fully discrete scheme. Now we describe the fully discrete scheme based on the LI 
approximation (1.3): find Uh € Xh for n = 1, 2,..., N 

n—1 

(2.8) (bol + T°‘Ah)UZ = - bj)U]f-^ + r“K, 

with Uh = Vh and = Phf{tn)- We focus on the homogeneous case, i.e., / = 0, and use 
the technique in [17] and [18] for convergence analysis. Throughout, we denote by 

OO 

^(0 = 

1=0 

the generating function of a sequence To analyze the fully discrete scheme (2.8), 

we first derive a discrete analogue of the solution representation (2.7). The fully discrete 
solution Wh := Ufi — Uh satisfies the following time-stepping scheme for n = 1, 2,..., N 

L^iWh) + AhWj: = -AhVh, 

with Wh = 0. Next multiplying both sides of the equation by and summing from 1 to 
OO yields 

°° ~ r 

Y. L'liWh)C + AhWhiV) = —rYp^hVh. 

Tl=l ^ 

Now we focus on the term Ti'(lF/,)5". By the definition of the difference operator 

I/", we have 

OO OO / n—1 \ 

^ L-,iWh)e = ^■“E + E(^i -c 

n=l n=l \ J = 1 / 

OO / n — 1 \ OO / n — 1 \ 

= E EU" -E E e 

ra=l \1 = 0 / ra=l \1 = 1 / 

:= 7 - II. 

Using the fact Wh = 0 and the convolution rule of generating functions (discrete Laplace 
transform), the first term 7 can be written as 

OO 71 

I = E(E = r-mwhio- 

71=1 J = 0 
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Similarly, the second term II can be written as 

oo n oo n—1 

II = r-“ 

n=l j = l n=l j = 0 

Hence, we arrive at 


= r-“(l - C)6(C)VF4C). 


Next we derive a proper representation for 6(^): 

1 " 


K0 = 




1-c 




cr(2 - a) ^ 


(i-OLu-i(e) 

cr(2 - a) ’ 


where Lip( 2 :) denotes the polylogarithm function defined by (see [15]) 

OO j 

Lip(«) = E 

The polylogarithm function Lip( 2 ) is well defined for \z\ < 1, and it can be analytically 
continued to the split complex plane C \ [l,oo); see [6]. With 2 = 1, it recovers the 
Riemann zeta function C,{p) = Lip(l). Therefore, the fully discrete solution Wh{^) can be 
represented by 

n2 \ -1 


Wh{^) = - 


€ 


(1-0= 


1 - ^ VCi'“r(2 - a) 


LiQ-i(0 + 2l/i AhVh- 


Simple calculation shows that the function Wh{C) is analytic at ^ = 0. Hence the Cauchy 
theorem implies that for q small enough, there holds 

\2 \ -1 


W^ = - 


-f 

2711 Jljl 


(1-cr 


27rii|5l^^ (1-OU VC^“r(2-a) 
Upon changing variable ^ = e~^'^, we obtain 


Lia-i(0 + ^ii AhVhd^- 


wr = —. 


' 

27ri Jpo 


1 — e“ 


-Lie —1 (e ) + 


AhVh dz, 


g-2Tpar(2 - a) 

where the contour T'^ := {z = — ln(£i)/r+ ii/ : jyj < vr/r} is oriented counterclockwise. By 
deforming the contour T® to T^ := {z G : |3(2:)| < ir/r} and using the periodicity of 
the exponential function, we obtain the following alternative representation for Wf^ 


/p^ 1 - Ve“"'^i'“r (2 - a) 

This representation forms the basis of the error analysis. 


(2.9) W." = 


LiQ-i(e ^^) + Hh AhVhdz. 


3. Error analysis of the fully discrete scheme 

In this part, we derive optimal error estimates for the fully discrete scheme (2.8). The 
analysis is based on the representations of the semidiscrete and fully discrete solutions, 
i.e., (2.7) and (2.9). Upon subtracting them, we may write the difference between lU^ 
and Whitn) as 

Wh{tn)-w;: = i+ii, 
where the terms I and II are given by 

I = f Ki{z)vhdz 

■'re,j\r^ 
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and 

//= 7 ^ [ K2{z))vhdz. 

2711 

Here the kernel functions Ki (z) and K 2 (z) are given by 

(3.1) ^i(z) ■= -z~^(z‘^ + Ah)~^Ah =■■ z~^Bi(z) 

and 

(3-2) J^2(z) := -^3^317 + A,^ ' A, =: 

respectively, and the auxiliary function i/j is define by 

^ r{2-a) 

Since the function is uniformly bounded on the contour Ft, we deduce 

||7Fi(^) - e-^-7F2(^)|| < \e-^l\\K^iz) - K2iz)\\ + |1 - e-^nil^i(^)ll 
(3.3) < c\\Ki(z) - K 2 iz)\\ + c\z\r\\Ki{z)\\ 

< c\\Ki{z) - K 2 iz)\\ +cr, 

where the last line follows from the inequality, in view of (2.3): 

||7Fi(.)l| = - I + z-iz- + A,)-^\\ < c|^r\ 

Thus it suffices to establish a bound on ||7Fi(2:) — 7^2(2)11, which will be carried out below. 
First we give bounds on the function x{^) = 

Lemma 3.1. Let x{z) =r“^(l —e“^’^). Then for all z € F^, there hold for some ci, C 2 > 0 
\x{z) - z\ < c\zf T and ci\z\ < \x(z)\ < C 2 \z\. 


Proof. We note that \z\t < tt/ sin0 for 2 : € Ft. Then the first assertion follows by 


|X(2) -2| < \zfr 


E 

j=0 


(i + 2)! 


< clzl r for 2 G Ft. 


Next we consider the second claim. The upper bound on x{^) is trivial. Thus it suffices 
to verify the lower bound. To this end, we split the contour Ft into three disjoint parts 
Ft = Ft U Ft U F^ , with Ft and FF being the rays in the upper and lower half planes, 
respectively, and F): being the circular arc. Here we set ^ = —zt with p = |^| € (0, tt/ sin 9). 
We first consider the case of 2 G Fjl, for which ^ = pe p G (l,7r/ sin^). Using the 

trivial inequality |cos(psin0)| < 1, we obtain 


(3.4) 


1 - e-"" 


— 1 

|e-^“"®cos(psin6l)-l- 

ig-pco.esin(psin6l)| 

ZT 



P 


/ -2pcos0 I 1 _ 9 -pcos6l\l/2 -pcosS _ i 

> ^ -+_i_^e- }_ ^ e -i ^ 


P P 

due to positivity and monotonicity of — l)/p as a function of p over the interval 

(0, 00 ). The case of 2 G FF follows analogously. Last, we consider 2 G Ft, the circular 
arc. In this case, by means of Taylor expansion, we have 


X{z) = z ( 1 + ^(-1)^ 

J = 1 


U + 1 )! / ' 


From this and the fact that p = | 2 r| < 1, it follows directly that |x(z)l > c\z\ for 2 G F)!. 
This completes the proof of the lemma. □ 
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Then by the trivial inequality ||5i(2:)|| < c and Lemma 3.1, we deduce that 
\\K,{z) - K,{z)\\ < Iz-^ - x(^)-'|||Bi(^)|| + \x{z)\-^\\Bi{z) - B 2 {z)\\ 

(3.5) < +c\z\-^\\Bi{z)-B 2 iz)\\ 

<CT + c\z\~^\\Bi{z) - B 2 {z)\\. 

Thus it suffices to establish a bound on ||_Bi( 2 ) — 52 ( 2 )||. This will be done using a series 
of lemmas. To this end, first we recall an important singular expansion of the function 
Lip(e“^) [6, Theorem 1]. 

Lemma 3.2. For p 7 ^ 1,2,.. the function Lip(e~^) satisfies the singular expansion 

°° , 

(3.6) Lip(e ^)~r(l-p) 2 *’ ^+ 5Z(“1) C(p-fc)-;g- as z ^ 0, 

k=0 


where C, is the Riemann zeta function. 

Remark 3.1. The singular expansion in Lemma 3.2 is stated only for 2 —>■ 0. However, 
we note that the expansion is valid in the sector T.e,5 [6, pp. 377, proof of Lemma 2]. 

We shall need the following results. The first result gives the absolute convergence of 
a special series involving the Riemann zeta function C). 

Lemma 3.3. Let \z\ < n/sinO with 9 £ (7r/2, 57r/6), and p = a — 1. Then the series 

(3.6) converges absolutely. 


Proof. Using the following well-known functional equation for the Riemann zeta function 
(see e.g., [14] for a short proof): for 2 ^ Z, there holds 

we obtain for p = a — 1 € (—1, 0) 

C(p-fc) =C{l-{l-p + k)) 

= (^2Tvy-p+k cos p + r{l-p + k)ai-p + k). 

By Stirling’s formula for the Gamma function r(a;), a; —>■ 00 [1, pp. 257] 
r(a; -I- 1) = (l -I- 0(a;“^)) 


and that <^(1 —p-ffc)—^-lasfc—cxa, we have 


lim 

k —>-oo 


k/KiP-k)\\z\'^ ^ 1 


V| 2 | < tt/ sin 6. 


k\ 2 sin 9 

Since for 9 G (7r/2, 57r/6), 2sin0 > 1, the series converges absolutely. 


□ 


Next we state an error estimate for the function 


-iP{zt) with respect to 2 “. 


Lemma 3.4. Let ip{z) = j^^^^hia-i{e ^). Then for the choice 9 £ (7r/2, 57r/6), there 
holds 

,1-e-"" 


-iP(zt) — Z°‘\ < c\z\ T 


Mz £ T.r. 
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Proof. Upon noting the fact 0 < \zr\ < tv/ sin^ and using Taylor expansion and (3.6), we 
deduce that for 2 G Fr, there holds 


{ztY 




Li,_i(e-"") = r(2 - a)(zr)“-" + ^(-1)^(1 - a - k) 


fc=0 


{jzt 

k\ 


Hence, the function can be represented by 


= E ^ 


j=i 




i-lfa-a-k) {zrf 


r(2-a) 


k\ 


= E 


{zr) 


a+j-2 oo 


{zry (—a — k) {zrY 


■ + EEi^E 




^ j! ^ r(2-a) 

= 1 ■' fc =0 ^ ' 


k\ 


and 


-y{zt) = T 

1=0 


{1 + 1)1 


^ (£E g (-l)"C(-« - k) {zr)'‘ 


J=i 


?! r(2 — a) k\ 

j=l k=0 ^ ' 


ct , C(<^ ~ 1) 2 2 -q , /-V/ 2+a 2\ 


r (2 - a) 

In view of the choice Q £ (7r/2,57r/6) and Lemma 3.3, the bound is uniform, since the 
series converges uniformly for 2 : € r^. Consequently, 


I-!—£- iP{zt)- z°‘\<\zYr'^ 


C(c^-l) 
r (2 - a) 


+ 0 {{zr)‘^)] <clzYr^-r 


from which the desired assertion follows. □ 

The next result gives a uniform lower bound on the function tp{z) on the contour r,-. 

Lemma 3.5. Let '^ken for any 9 close to '7r/2, there holds for 

any 5 < 7r/2r 

l'*/’( 2 r)| > c > 0 G Ft. 

Proof. Since for z G F.^, P-^l < t^/t and 2 ^ ( — 00 , 0], by [18, Lemma 1], there holds 


Yiz) =Cc [ 
Jo 


1 - e-^-‘ 


■ ds, 


with the constant Ca = sin(7r(l — a))/7r. To prove the assertion, we again split the contour 
Ft into F = Ft U Ft U Ft , where Ft and F 7 are the rays in the upper and lower half 
planes, respectively, and F^ is the circular arc of the contour Ft, and discuss the three 
cases separately. We first consider the case 2 G F^f and set zt — pe'® = pcos9 + ip sin 0 
with S < p < Tr/sinS. Upon letting r = pcos9 and rf = psinS then 

„a — l 


%l){zr) = Ca 

Jo 

r 

= Ca 

Jo 


1 — e 


1 - ^-(r++)-a g 


■ ds 


■ ds 


I 


1 — e * cos <j> + ie ’’ sin 
°° — e“'’)(l — e“’'“'’cos0 — ie“’'“'’sim^) 


(1 — 6“’'“® cos 0)2 + e' 
It suffices to show that the real part 


ds. 


lfiip{zr) = 


)=cc. r 

Jo 


ds 


1 — 2e“’ ^ cos 0 + e“2’ 2 s 
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is bounded from below by some positive constant c. First we consider the case <j> = p sin 6 € 
[ 71 / 2 , tt], for which cos0 < 0 and thus 

0 < 1 — cos < 1 — cos <j) 

^1 o —r — s I I —2r — 2s 

< 1 — 2e cos 0 + e 


Consequently, 

POO 

>Ca — e”'*) ds = co- 

Jo 

Next we consider the case (f> € (0, tv/2), for which costf) > 0. Further we fix 0 = n/2, and 
thus r = pcosO = 0 and cos cf) = cos(psin0) = cosp > 0. Then 

1 — e~' cos (j) > 1 — e~° and 0 < 1 — 2e~’ ® cos 4> + e~^' < 2, 


and accordingly the real part Ki/;(zr) simplifies to 



s“ ^(1 —e ®)(1 —e “cosp) 
1 — 26“® cos p + 

— e~“)^ds > Cl. 


ds 


Then by continuity of iRi/)( 2 r), we may choose an angle 6 € (ir/S, Stt/G) such that for any 
z £ r^, there holds K'(/>(zt) > C 2 . Repeating the above argument shows also the assertion 
for the case z G F^. It remains to show the case z G r“. For any fixed p G (0, tt/2) and 
9 G [— 71 / 2 , 71 / 2 ], cos(j) = cos(psin6') > 0, r = pcosO > 0. Consequently 


1 — e “ cos 0 > 1 — e “ cos (/ > 1 — e 


and 

1 o —r — s I I —2r —2s ^ i , —2r —2s ^ r» 

1 — 2e cos0 + e <l + e <2. 


These two inequalities directly imply 

""-^(l-e-“)(l-e-"-“cos (j}) 


= c, 

> 


Jo ~ 


— 2e“'’“® cosd) + e 


-2r — 2s 


ds 


Cq: / a —2/1 —s\2 1 \ 

Y J s (1 - e ) ds>C3. 


Then by continuity, we may choose an angle 0 > 7r/2 such that for z G F/, there holds 

%iP{zt) > C4 > 0. □ 


The next result shows a “sector-preserving” property of the mapping xi (z)’- there exists 
some ^0 < 71, such that Xi{^) ^ for all 2 G Eg. This property plays a fundamental 
role in the error analysis below. 

Lemma 3.6. Let il>{z) = p'[^~^) Lia-i(e~^) and X^i^) = — 'iji^zr). Then there exists 

some 6o G ( 71 / 2 , ti) smc/i that Xi(^) € Eg^ for all z G Eg. 


Proof. Like before, for zt = pe‘®, we denote by r = pcos9, cj> = psin0 and Ca = sin(7r(l — 
a))/7r. Then the real part Kxi(^) ^^ncl the imaginary part £yxi('2) of the kernel Xi(^) S'l'e 
given by 


(3.7) 

and 


KXi(2) 



s°-2(i _ e-^)(i 6-^”-“ - 6-”-“ cos (t> - e"” cos (f) 

1 — 26“’"“^ cos <j) + 


ds 


(3.8) 


®Xi(«) 



s°-^(l-e-“)^e-”sin0 
1 — 2e“’ “ cos <j) + ’ 
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respectively. Obviously, for 6 € [—7r/2,7r/2] then r = pcosO >0, 0<e ^^<1 and thus 

H I —2r — s —r — s , —r j. \ , —2r — s —r — s —r 

1 + e — e cos (p — e cos 0>l + e —e —e 

s 0 < 1 - 2e' 

Cc(l-e-") 


Meanwhile, r > 0 implies 0 < 1 — 2e ^ ” cos </> + e < 4, and consequently 




4r“ 


/■ 


\l-e-rds=^{l-e-n, 


with the constant ^ — e~‘‘)^ds. 

Next we consider the case \8\ > -k jl. It suffices to consider the case 8 > 7r/2, and the 
other case 8 < —-k 12 can be treated analogously. Let zt = pe‘® with p € (0, oo). First, 
clearly, for ij> = psind € [7r/2,7r], cos0 < 0, and there holds 


0 < 1 - 26-’'-'’ cos 0 + e-"’'-"'’ < (1 + e-’'-'’)". 


and thus 




r“(l + e-'’)' 


POO 

■ / s“-^(l — 6-'’) ds > 0. 

Jo 


Second, we consider the case </> = psind € (0,7r/2). There are two possible situations: (a) 
1 — e~'" cos (^ > 0 and (b) 1 — e~'" cos (j) < 0- In case (a), we have 

1 + e — e cos (p — e cos (p 
>1 + cos^ (p — cos(p — e~^ costp 

=(1 — 6-’^-'* cos </>)(! — e~'" coscp) > (1 — e-'’)(l — e~'" coscp). 


Consequently, 


Kyi ( 2 ) > 


Ca(l — e ’’ cos (p) 


POO 

Jo ~ 


'( 1 - 6-^)2 


■ ds > 0. 


— 2e“^“® cos (p + 

In case (b), we may further assume 5fixi(2) < 0, otherwise the statement follows directly. 
Then appealing to (3.7) and using the trivial inequality |cos())| < 1, we deduce that 
e~'~ > 1 and 

1 — e~^ cos (p < 0 and — e~'" cos (p > 0. 

With the help of these two inequalities, and the assumption ( 2 ) < 0, we arrive at 
0 > 1 + cos cp — e~^ cos cp 

= (1 — e~^ cos (p) + e-“(e-^’' — e~^ cos (p) 

> {cos(p — e~’") + — e~^ costp) 

> 6-’’(cos cp — e~’") + — e~’" cos p) 

= (1 — cosp — 6-^’’), 

where the first and fourth inequalities follow from SRxi(z) < 0 and > 1, respectively. 
Consequently, 


|5Hxi( 2)| < -^(e ^’■-6 ^cosp), 


with the constant 


Cl = Ci(r, p) = 


, P)= Ca 

Jo 


'(l-e-“)2 


1 — 2e * cos p + e 
Meanwhile, it follows directly from (3.8) that 


■ ds. 


Pxi(2)| = -^e ’’sin.; 


Therefore, 


|gIXi(2)| ^ sin(psin6>) 

|5Rxi(2)| ~ 6“^“®® — cos(psind) 


=: gip)- 
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Now set gi{p) = sin(psin 6 ) and 32 (p) = e —cos(psinP). Since for p G (0, 7 r /(2 sin0)) 
and 6 > 7r/2, 

linip(p) =-tan 6 l, pi(p),P 2 (p) > 0 , pj(p) < 0 and P 2 (p) > 0 , 

p^o 

i.e., the function g[p) is monotonically decreasing on the interval [ 0 , 7 r /2 sin 0 ], we deduce 
rn 'JJo ■ = >0. 

pE(0,7r/(2 sin 0)) 

This completes the proof of the lemma. □ 


Remark 3.2. By the proof of Lemma 3.6, the sector Egp depends on the choice of the 
angle 9. For 9 —>■ tt/ 2, it is contained in the sector for any e > 0. 

Lemma 3.7. Let 9 be close to 7r/2, and 5 < 7r/2r. Then for Ki(z) and K 2 {z) defined in 
(3.1) and (3.2), respectively, there holds 

\\Kl{z) - K 2 {z)\\ < CT \Jz&Vr. 

Proof. For the operators Bi{z) and B 2 {z) defined in (3.1) and (3.2), respectively, we have 
Bi( 2 ) =-/+ 2“(z“ + Ah)“^ and B 2 {z) =-I + xi{z)ixi{^) + ^h)~^, 
where Xi('*) = — if{zT). Then by Lemma 3.4 

\\B,iz) - B 2 {z)\\ < \xi{z) - 2“|||(xiW + + kriKxiW + An)-^ - (^“ + Ah)-^\\ 

< c\z\\^-uxi(z)+ u)-^ + knxi(^) - ^“iii(xi(^)+ A^r^z-+A^r^w 

<c\z\\^-\\{Xi{z) + Ah)-^\\. 

Now we note that 

\Xliz)\ = |x(z)|t^"“|V'(2t)| > c|2|t^"“. 

By Lemma 3.6, x(2^) € Sgg for some 9o G (7r/2 , tt). Thus by Lemma 3.5 and the resolvent 
estimate (2.3), we have 

||A'i( 2 ) — ^^2(2)11 < c\z\~^\z\'^T^~°‘\xi(z)\~^ +CT < CT, 
and the desired estimate follows immediately. □ 


Now we can state an error estimate for the discretization error in time for nonsmooth 
initial data, i.e., v G L^{Q). 

Theorem 3.1. LetUh and he the solutions of problems (2.5) and (2.8) withv G L^(fl), 
Uh = '^h ~ PhV and / = 0, respectively. Then there holds 

||'t^/i(tTi) Uf, ||i:.2(n) — 


Proof. It sumces to bound the terms I and 
arrive at the following bound for the term I 


•TTf {t sin 6) 


rtji cos 6 


dr + 




ll-^lli,2(n) < cr\\vh\\L‘i{n) / 

(3.9) 

< ct"V||n/i||4,2(n). 

By Lemma 3.7 and direct calculation, we bound the term II by 

roo 

l7r/(TsinS) 

< cr\\vh\\L'^(n) f dr < crt;(^||n/i||4,2(n). 

Jo 

Combining estimates (3.9) and (3.10) yields \\wh{tn) — 14^11^2(0) < CTtn^\\vh\\L'^(n) and 
the desired result follows from the identity Lffi — Uh{tn) = WJf — Wh{tn) and the stability 
of the projection Ph in L/^iQ,). □ 
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Remark 3.3. The L^{Q,) stability of the LI scheme follows directly from Theorem 3.1. 

Next we turn to smooth initial data, i.e., v € D(A) = H^{Q) nRo(f2). To this end, we 
first state an alternative estimate on the solution kernels. 

Lemma 3.8. Let 6 be close to tv/ 2, and S < 'k/2t. Further, let Kl{z) = +Ah)~^ 

and K 2 {z) = + Ah)~^. Then for any z €'Ss,e, there hold 

\\Kl{z)-K^ 2 {^)\\<c\z\-r. 

Proof. Like before, we set Bf{z) = —( 2 “ + Ah)~^ and B^iz) = —{xi{^) + Ah)~^. Then 
by Lemmas 3.4 and 3.5, we have 

\\Bt{z) - B^iz)\\ < \xiiz) - 2 “|||( 2 “ + A,)-i||(xi( 2 ) + A,)-^\\ 
<c| 2 |V-“| 2 r“|xi( 2 )ri <c|z|i-“r. 

The rest follows analogously to the derivation (3.3). □ 


Now we can state an error estimate for Av € L^{Q). 

Theorem 3.2. LetUh and UJ^ be the solutions of problems (2.5) and (2.8) withv € 

Uh = I’h ~ RhV and / = 0, respectively. Then there holds 

\\uh{tr^) - Uh\\i^2(^n) < crt“"^||Tn||^2(n). 

Proof Let Kt{z) = -z-'^{z°‘ + Ah)-^ and Ki{z) = -xiz)-\xi{z) + Ah)-\ Then we 
can rewrite the error as 

e’'^” Kl{z)AhVhdz 

\rx 

e^*”(7i'i(2) - Kl{z))AhVhdz = 7 + 7J. 

By Lemma 3.8 we have for z € Ft 


(3.11) 


Wu{u)-w;:~ f 

27ri ./ 


||7f?(2)-e-^"7^2=(2)||<c|2r“r. 

By setting S = Iftn and for all z £ we derive the following bound for the term 77 


P^lln2(n) < CT\\AhVh\\L2(ei) 


pTTI {t sin 6 ) 


(3.12) 

< ctZ~^r\\AhVh\\L^(n)- 
Now (2.3) implies that for all z (zFs.e 


rtn cos 9 — 

e r 




cosV^.a-l . , 


li,2(n) < c\\AhVh\\L'^{xt) 


(3.13) 


/' 


rtrj cos 0 —a — 1 7 

e r dr 


/(r sin 9) 

POO 

< CT\\AhVh\\L2(n) / dr < crt““^||T,in^||^2(n). 

Jo 


Then the desired result follows directly from (3.12), (3.13) and the identities Uf — Uh{tn) = 
W" - Wh{t„) and AhRh = PhA. □ 


Remark 3.4. The convergence behavior of the LI scheme is identical with that for the 
convolution quadrature generated by the backward Euler method, which also converges at 
an 0{r) rate [11]. In particular for smooth initial data v G L){A), the time discretization 
error by both schemes contains a singularity This singularity reflects the limited 

smoothing property of the solution u [20, Theorem 2.1] 

l|Sr'“(i)IL2(n) < c[[An[[^2(n), 

whereas the first order derivative u'{t) is unbounded at t — 0. 
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Example 3.1. To illustrate the convergence rate in Theorem 3.2, we give a trivial ex¬ 
ample. Consider the following initial value problem for the fractional ordinary differential 
eguation: 

df u + u = Q, Vi > 0, with u(0) — 1. 

The exact solution u att = r is given by u{t) = Ea,i{—T°‘), where Ea,i(z} = 2 *^/r(Q:fc+ 

1) is the Mittag-Leffier function. For small t, the LI scheme at the first step is given by 

OO 

= (1 + r(2 - a)r°‘)~^ = 1 + ^(-l)"(r(2 - Q)r“)’". 

n=l 

Then the difference between and u{t) is given by 

u{t) — U^ = (r (2 — a) — r(a + + Ctt'^°‘, 

with Ct = “ r(2 — Since \ct\ < co for small t, we 

deduce that 

|M(r) — (7^1 ~ = tf~^T. 

This confirms the convergence order in Theorem 3.2. 

Last, the error estimates for the fully discrete scheme (2.8) follow from Theorems 2.1, 

3.1 and 3.2 and the triangle inequality. 

Theorem 3.3. Let u and Ufi be the solutions of problems (1.1) and (2.8) with = Vh 
and f = 0, respectively. Then the following estimates hold. 

(a) If V £ D{A) and Vh = RhV, then for n > 1 

\\uit„) - 17^111,2(0) < + h'^)\\Av\\L2^ci). 

(b) If V £ L^{Ll) and Vh = PhV, then for n > 1 

||u(in) — < c{rt.„^ + li^iTi°')ll'y||i,2(r;)- 

Remark 3.5. For v £ D{A), we can also choose Vh = PhV by the stability of the LI 
scheme. Hence, by interpolation we deduce 

Mtn) - < c(ri7'+“" + 0 < a < 1. 

4. Time-space fractional differential problem 

The convergence theory developed in Section 3 may be extended to more general secto¬ 
rial operators A, i.e., (a) The resolvent set p{A) contains the sector {z : 6 < \ arg 2 | < tt} 
for some 6 £ (0,7r/4); (b) ||(2;7 -|-7l)“^|| < M/\z\ for 2 £ T,Ty-g and some constant M. 
The technical restriction 9 £ (0,7r/4) stems from Remark 3.2. This in particular covers 
the Riemann-Liouville fractional derivative of order /3 € (3/2,2); see Lemma 4.1 below. 
Specifically, we consider the following one-dimensional space-time fractional differential 
equation 

dtU-^D^u = f, inn = (0,1) T>t>0, 

(4.1) u = 0, on dLl T > t > 0, 

u(0) = V, in LI, 

with OL £ (0,1) and /3 G (3/2, 2). Here with n— l</3<n, n€N, denotes the 
left-sided Riemann-Liouville fractional derivative qD^u of order j3 defined by [13, pp. 70]: 

(4.2) ix-sr-^-Ms)ds. 

The right-sided version of Riemann-Liouville fractional derivative is defined analogously 
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The model (4.1) is often adopted to describe anomalous diffusion process involving both 
long range interactions and history mechanism. 

The variational formulation of (4.1) is to find u € such that (see 

[8]) 

(4.3) {dTu, ip) + Aiu, p) = if, p) yp G 
with w(0) = V, where the sesquilinear form A(-, ■) is given by 

A{p,ip) = - (oD^^^p, . 

It is known (see [.5, Lemma 3.1] and [8, Lemma 4.2]) that the sesquilinear form ^(-j •) is 
coercive and bounded on the space Then Riesz representation theorem implies 

that there exists a unique bounded linear operator A : such that 

Define D{A) = {i/i G : Aip G L^{Q)} and an operator A : D{A) —>■ L^{Q) by 

(4.4) Aip, V>) = {Ap, ^), pe D{A), ^ G 

We recall that the domain D{A) consists of functions of the form off / — (off/)(1)®^~^, 
where / G L^{Q,) [8]. The term G fff“^"'"^(fl), <5 G [1 — /3/2,1/2), appears because 
it is in the kernel of the operator . The presence of the term indicates that the 
solution u usually can only have limited regularity. 

Lemma 4.1. For j3 G (3/2,2), the resolvent set p{A) contains the sector {z : 9 < 

I a'rg( 2 )| < tt} for any 6 £ {{2 — /3)7r/2,7r/4). 

Proof. Let u be the zero extension of u. Recall that for s G (1/2, 1) and Q, = (0,1), 

ll«ll5nn) = IIICr-^(“)(0llL2(K), 

where F{u) is the Fourier transform of m, is a consistent and well-defined norm on ff“(Sf). 
Further, we note that for u G ff^^^(fl) (see [.5] and [8]) 

3f(A(M, u)) > co||u||^,3/2(n) with co = cos((2 - /9/2)7r). 

Further, we recall the fact that for u G Co°(fl) and s G (1/2,1) 

\\oDxU\\]^2(^q') = II -off/cW||i;,2(Q) < II -ooffxW||i2(j{) 

= l|-f( -offffw)||^2(]g) = ||f^( -off//u)||l,2(]g) 

= iiier-F(«)(oiiL2(«) = 11^^115.(0)- 

By the density of C^{Q) in for s G (1/2, 1) we obtain for all u G ff^'^^(fl) 

llof/f^ u||£2(o) < ||m||^3/2(o). 

Likewise, the right sided case follows: 

'w||L2(n) ^ Il'*^ll53/2(n)- 

Thus for u,v £ ff^^^(D), there holds 

|24(m,u)| < Cl||u||^4i/2(n)||u||ifAi/2(n) with Cl = 1. 

Then by Lemma 2.1 of [9], we conclude that the resolvent set p{A) contains the sector 
{z : 9 < \ argzl < tt} for all 9 £ ((2 — I3I2)'x,'k 12). In particular, for fi > 3/2, we may 
choose 0 G ((2 —/3/2)7r, 7r/4). □ 
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By Lemma 4.1 and Remark 3.2, we can apply the theory in Sections 2 and 3 to derive 
a fully discrete scheme based on the LI scheme in time and the Galerkin finite element 
approximation in space. First we partition the unit interval into a uniform mesh with 
a mesh size h = 1/M. We then define Vh to be the set of continuous functions in V which 
are linear when restricted to the subintervals [xi,Xi+i], i = 0,... ,M — 1. Further, we 
dehne the discrete operator Ah ’■ Vh ^ Vh by 

{Ah^p, x) = A{ifi, x) V(p, X G Vh. 

The corresponding Ritz projection Rh : Vh is dehned by 

(4.5) A{Rh^P, x) = Ai^, x) x&Vh. 

Then the fully discrete scheme for problem (4.1) based on the LI approximation (1.3) 
reads: hnd 11// for n = 1, 2 ,..., 

n 

(4.6) (7 + r^AhW// = b^U°h + ^(5,-1 - h,)Ur^ + 
with Uh = Vh and F// = Phf{t„). 

Last, we state the error estimates for the fully discrete scheme (4.6). These estimates 
follow from Theorems 3.1 and 3.2 and the error estimates for the semidiscrete Galerkin 
scheme, which can be proved using the operator trick in [9], and thus the proof is omitted. 

Theorem 4.1. Let u and UJ/ be the solutions of problems (4.1) and (4.6) with Uh = Vh 
and / = 0, respectively. Then for 5 € [1 — P/2, 1/2) the following estimates hold. 

(a) If V £ D{A) and Vh = RhV, then for n > 1 

l|M(fn) - t/;/||i2(n) < c(rt““^ + h^~^~^^)\\Av\\h2^ny 

(b) If V £ L^{Q) and Vh = PhV, then for n > 1 

Wu/tr,) - [/^||L2(n) < c{Tt~^ + h^~^*^t~°‘)\\v\\h2^n)- 

5. Numerical experiments and discussions 

Now we present numerical results to verify the convergence theory in Sections 3 and 
4, i.e., the 0{t) convergence rate. We consider the subdiffusion case (1.1) and time-space 
fractional case (4.1) separately. For each model, we consider the following two initial data: 

(a) = (0, 1), and v = sin(27ra;) € 77^(12) n 77o(f2). 

(b) 72 = (0,1), and v = £ 77^/'*“'(72), with e £ (0,1/4); 

We measure the error e" = u{tn) — Uh by the normalized errors ||e"||/^ 2 (n)/||u||i^ 2 (f;). In the 
computations, we divide the unit interval 72 = (0,1) into M equally spaced subintervals 
with a mesh size h = 1/M. Likewise, we fix the time step size r at r = t/N, where t is 
the time of interest. In this work, we only examine the temporal convergence rate, since 
the space convergence rate has been examined earlier in [10, 9]. To this end, we take a 
small mesh size h = 2~^^, so that the spatial discretization error is negligible. 

5.1. Subdiffusion. The exact solution can be written explicitly as an infinite series in¬ 
volving the Mittag-Leffler function Ea,p{z) dehned by 

oo h 

see [10] for the details. The Mittag-Leffler function Ea,p{z) can be evaluated efflciently 
by the algorithm developed in [21]. The numerical results for cases (a) and (b) are shown 
in Table 2, where rate refers to the empirical convergence rate, and the number in the 
bracket refers to the theoretical rate from Theorem 3.3. The results fully conhrm the 
theoretical prediction: for both smooth and nonsmooth data, the fully discrete solution 
Uh converges at a rate 0(r), and the rate is independent of the fractional order a. Further, 
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for fixed t, the error increases with the fractional order a. This might be attributed to 
the local decay behavior of the solution u: the larger is the fractional order a, the slower 
is the solution decay around t = 0. According to Remark 3.5, we have 

HU) - [/."IIlrh) < ^ e [o, i]. 

Thus the temporal error deteriorates as the time tn —>■ 0 at a rate like and for 

Av G L^(Q.) and v G D{A^^^~'^), with e G (0,1/8), respectively. In particular, for fixed 
N, the error behaves like and for cases (a) and (b) with a = 0.5, respectively. 
This is clearly observed in Table 3, thereby showing the sharpness of the error estimates. 
Further, we observe that the LI scheme fails to converge uniformly (in time) at a first 
order even though the initial data v in case (a) is very smooth, i.e., v G D{A^) for any 
p > 0, cf. Table 3. This numerically confirms the observation in Remark 3.4. 


Table 2. The L^-norm of the error for the subdiffusion equation with 
initial data (a) and (b) at t = 0.1 with h = 2“^^, r = 1/N. 


a 

N 

10 

20 

40 

80 

160 

320 

rate 

0.1 

(a) 

1.46e-4 

7.18e-5 

3.55e-5 

1.77e-5 

8.82e-6 

4.40e-6 

« 1.01 (1.00) 


(b) 

3.95e-4 

1.93e-4 

9.57e-5 

4.76e-5 

2.38e-5 

1.19e-5 

« 1.00 (1.00) 

0.5 

(a) 

1.22e-3 

5.89e-4 

2.88e-4 

1.43e-4 

7.08e-5 

3.52e-5 

« 1.01 (1.00) 


(b) 

3.65e-3 

1.73e-3 

8.36e-4 

4.09e-4 

2.02e-4 

l.OOe-4 

« 1.02 (1.00) 

0.9 

(a) 

7.01e-3 

3.05e-3 

1.39e-3 

6.53e-4 

3.12e-4 

1.50e-4 

« 1.07 (1.00) 


(b) 

1.54e-2 

7.67e-3 

3.79e-3 

1.87e-3 

9.23e-4 

4.55e-4 

« 1.02 (1.00) 


Table 3. The L^-norm of the error for the subdiffusion equation with 
initial data (a) and (b) with a = 0.5, h = 2“^® and N = 10. 


t 

le-5 

le-6 

le-7 

le-8 

le-9 

le-10 

rate 

(a) 

2.94e-3 

1.05e-3 

3.45e-4 

l.lle-4 

3.51e-5 

1.1 le-5 

« 0.50 (0.50) 

(b) 

3.02e-3 

2.56e-3 

2.18e-3 

1.86e-3 

1.58e-3 

1.35e-3 

« 0.07 (0.06) 


5.2. Time-space fractional problem. Now we present numerical results for the space 
time fractional problem. Since the exact solution is not available in closed form, we 
compute the reference solution by the second-order backward difference scheme from [11] 
on a much finer mesh, i.e., N = 1000. The numerical results for case (a) with different 
a and j3 values are presented in Table 4. A first-order convergence is observed, and the 
convergence rate is independent of the time- and space-fractional orders. Interestingly, the 
observation is valid also for the case = 5/4, for which the theory in Section 4 does not 
cover, awaiting further study. For a fixed a value, the error decreases with the increase 
of the fractional order P, which indicates that the solution decays faster as P approaches 
two. This is also consistent with the fact that the closer is the P value to unity, the more 
singular is the solution, and thus more challenging to approximate numerically, cf. [8]. 
For the nonsmooth case (b), we are particularly interested in the case of small t. Thus 
we present the numerical results for t = QA, t — 0.01 and t = 0.001 in Table 5. The 
experiment shows the first order convergence is robust for nonsmooth data even if t is 
close to zero. Like before, for fixed n, let tn —>■ 0, the error behaves like for case (a), 
which is fully confirmed by the numerical results, cf. Table 6, indicating the sharpness of 
the estimate in Theorem 4.1. 
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Table 4. The L^-norm of the error for the time-space fractional prob¬ 
lem with initial data (a) at t = 0.1 with h = 2“^®, r ^ 0. 


a 

P\N 

5 

10 

20 

40 

80 

rate 


1.25 

1.37e-3 

6.55e-4 

3.21e-4 

1.59e-4 

7.90e-5 

« 1.01 (1.00) 

0.1 

1.5 

8.41e-4 

4.03e-4 

1.98e-4 

9.78e-5 

4.87e-5 

« 1.01 (1.00) 


1.75 

5.08e-4 

2.44e-4 

1.19e-4 

5.92e-5 

2.94e-5 

« 1.01 (1.00) 


1.25 

1.52e-2 

6.69e-3 

3.12e-3 

1.50e-3 

7.32e-4 

« 1.03 (1.00) 

0.5 

1.5 

8.22e-3 

3.70e-3 

1.75e-3 

8.49e-4 

4.17e-4 

« 1.05 (1.00) 


1.75 

4.69e-3 

2.14e-3 

1.02e-3 

4.97e-4 

2.45e-4 

« 1.03 (1.00) 


1.25 

6 .01e-2 

3.19e-2 

1.62e-2 

8.07e-3 

3.99e-4 

« 1.01 (1.00) 

0.9 

1.5 

5.61e-2 

2 .86e-2 

1.42e-2 

6.95e-3 

3.39e-3 

« 1.03 (1.00) 


1.75 

3.58e-2 

1 .66e-2 

7.74e-3 

3.66e-3 

1.75e-3 

« 1.07 (1.00) 


Table 5. The L^-norm of the error for the time-space fractional prob¬ 
lem with initial data (b) with /3 = 1.5, h = 2“^^, t = 0.1, 0.01 and 
0 . 001 . 


a 

t\N 

5 

10 

20 

40 

80 

rate 


0.1 

1.53e-3 

7.32e-4 

3.59e-4 

1.77e-4 

8.83e-5 

« 1.01 (1.00) 

0.1 

0.01 

1.74e-3 

8.34e-4 

4.08e-4 

2.02e-4 

l.Ole-4 

« 1.01 (1.00) 


0.001 

1.94e-3 

9.31e-4 

4.56e-4 

2.25e-4 

1.12e-4 

« 1.01 (1.00) 


0.1 

1.39e-2 

6.36e-3 

3.01e-3 

1.45e-3 

7.11e-4 

« 1.04 (1.00) 

0.5 

0.01 

1 .22e-2 

5.86e-3 

2.85e-3 

1.40e-3 

6.89e-4 

« 1.03 (1.00) 


0.001 

8.02e-3 

3.83e-3 

1.86e-3 

9.12e-4 

4.50e-4 

« 1.02 (1.00) 


0.1 

1.99e-2 

1 .02e-2 

5.15e-3 

2.60e-3 

1.31e-3 

« 1.00 (1.00) 

0.9 

0.01 

1 .21e-2 

6.10e-3 

3.05e-3 

1.52e-3 

7.53e-4 

« 1.00 (1.00) 


0.001 

7.63e-3 

3.83e-3 

1.91e-3 

9.51e-4 

4.73e-4 

« 1.00 (1.00) 


Table 6 . The L^-norm of the error for the time-space fractional prob¬ 
lem with initial data (a) and (b), a = 0.5 and /3 = 1.5, as t —>■ 0 with 
h = 2“^^ and N = 5. 


t 

le-5 

le-6 

le-7 

le-8 

le-9 

le-10 

rate 

(a) 

2.60e-3 

8.90e-4 

2.96e-4 

9.71e-5 

3.17e-5 

1.03e-5 

« 0.48 (0.50) 

(b) 

4.74e-3 

3.76e-3 

3.02e-3 

2.44e-3 

1.99e-3 

1.63e-3 

« 0.09 (—) 


6. Conclusions 

In this paper we have revisited the popular LI time-stepping scheme for discretizing the 
Caputo fractional derivative of order a € (0,1), arising in the modeling of subdiffusion, 
and rigorously established the first order convergence for both smooth and nonsmooth 
initial data. This result complements existing convergence analysis, which assumes a 
regularity in time. The extensive numerical experiments fully verify the sharpness of the 
estimates and robustness of the scheme. The convergence analysis is valid for more general 
sectorial operators, and in particular, it covers also the space-time fractional differential 
equations, for which we have also provided error estimates. 

In view of the solution singularity for time t close to zero, it is natural to consider the LI 
scheme on a nonuniform time mesh in order to arrive at a uniform first-order convergence. 
However, the generating function approach used in our analysis does not work directly 
in this case, and it is an interesting open question to rigorously establish error estimates 
directly in terms of the data regularity. A second interesting future research problem is 
the convergence rate analysis of the scheme for the diffusion wave equation, which involves 
a Caputo fractional derivative of order a € (1,2) in time. 
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